home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 10 / AACD 10.iso / AACD / Online / SpeakFreely / src / lpc / lpc.c next >
C/C++ Source or Header  |  2000-05-18  |  9KB  |  429 lines
  1. #include <stdio.h>
  2. #ifndef NO_LPC_FIX
  3. #include <stdlib.h>
  4. #endif
  5. #include <math.h>
  6. #include <sys/types.h>
  7. #include <netinet/in.h>
  8.  
  9. #include "lpc.h"
  10. #include "../ulaw2linear.h"
  11.  
  12. #if defined(Solaris) || defined(HEWLETT_PACKARD) || defined(__FreeBSD__)
  13. #define bcopy(a, b, n)      memmove(b, a, n)
  14. #endif
  15.  
  16. #define MAXWINDOW    1000    /* Max analysis window length */
  17. #define FS        8000.0    /* Sampling rate */
  18.  
  19. #define DOWN        5    /* Decimation for pitch analyzer */
  20. #define PITCHORDER    4    /* Model order for pitch analyzer */
  21. #define FC        600.0    /* Pitch analyzer filter cutoff */
  22. #define MINPIT        50.0    /* Minimum pitch */
  23. #define MAXPIT        300.0    /* Maximum pitch */
  24.  
  25. #define MINPER        (int)(FS/(DOWN*MAXPIT)+.5)    /* Minimum period */
  26. #define MAXPER        (int)(FS/(DOWN*MINPIT)+.5)    /* Maximum period */
  27.  
  28. #define WSCALE        1.5863    /* Energy loss due to windowing */
  29.  
  30. #define BUFLEN        ((FRAMESIZE * 3) / 2)
  31.  
  32. #define SILENCEFIX        /* Enable absolute silence fix */
  33.  
  34. /*  The LPC coder does truly awful things when driven into clipping.
  35.     If you set GAIN_ADJUST to a number less than 1.0, samples will
  36.     be scaled by that factor to avoid overdriving the coder.  */
  37.  
  38. #define GAIN_ADJUST 0.9
  39.  
  40. static struct lpcwork {
  41.     float w_s[MAXWINDOW], w_y[MAXWINDOW], w_h[MAXWINDOW], w_w[MAXWINDOW];
  42. } *work = NULL;
  43.  
  44. static float fa[6], u, u1, yp1, yp2;
  45.  
  46. #define lin_to_ulaw(x) audio_s2u(x)
  47. #define ulaw_to_lin(x) audio_u2s(x)
  48.  
  49. static void
  50. auto_correl(w, n, p, r)
  51.   float *w;
  52.   int n;
  53.   int p;
  54.   float *r;
  55. {
  56.     int i, k, nk;
  57.  
  58.     for (k = 0; k <= p; k++) {
  59.         nk = n - k;
  60.         r[k] = 0.0;
  61.         for (i = 0; i < nk; i++)
  62.             r[k] += w[i] * w[i + k];
  63.     }
  64. }
  65.  
  66. static void
  67. durbin(r, p, k, g)
  68.   float *r;
  69.   int p;
  70.   float *k;
  71.   float *g;
  72. {
  73.     int i, j;
  74.     float a[LPC_FILTORDER + 1], at[LPC_FILTORDER + 1], e;
  75.  
  76.     for (i = 0; i <= p; i++)
  77.         a[i] = at[i] = 0.0;
  78.  
  79.     e = r[0];
  80.     for (i = 1; i <= p; i++) {
  81.         k[i] = -r[i];
  82.         for (j = 1; j < i; j++) {
  83.             at[j] = a[j];
  84.             k[i] -= a[j] * r[i - j];
  85.         }
  86. #ifdef SILENCEFIX
  87. if (e == 0) {
  88.     *g = 0;
  89.     return;
  90. }
  91. #endif
  92.         k[i] /= e;
  93.         a[i] = k[i];
  94.         for (j = 1; j < i; j++)
  95.             a[j] = at[j] + k[i] * at[i - j];
  96.         e *= 1.0 - k[i] * k[i];
  97.     }
  98.  
  99. #ifdef SILENCEFIX
  100. if (e < 0) {
  101.     e = 0;
  102. }
  103. #endif
  104.     *g = sqrt(e);
  105. }
  106.  
  107. static void
  108. inverse_filter(w, k)
  109.   float *w;
  110.   float *k;
  111. {
  112.     int i, j;
  113.     float b[PITCHORDER + 1], bp[PITCHORDER + 1], f[PITCHORDER + 1];
  114.  
  115.     for (i = 0; i <= PITCHORDER; i++)
  116.         b[i] = f[i] = bp[i] = 0.0;
  117.  
  118.     for (i = 0; i < BUFLEN / DOWN; i++) {
  119.         f[0] = b[0] = w[i];
  120.         for (j = 1; j <= PITCHORDER; j++) {
  121.             f[j] = f[j - 1] + k[j] * bp[j - 1];
  122.             b[j] = k[j] * f[j - 1] + bp[j - 1];
  123.             bp[j - 1] = b[j - 1];
  124.         }
  125.         w[i] = f[PITCHORDER];
  126.     }
  127. }
  128.  
  129. static void
  130. calc_pitch(w, per)
  131.   float *w;
  132.   float *per;
  133. {
  134.     int i, j, rpos;
  135.     float d[MAXWINDOW / DOWN], k[PITCHORDER + 1], r[MAXPER + 1], g, rmax;
  136.     float rval, rm, rp;
  137.     float a, b, c, x, y;
  138.     static int vuv = 0;
  139.  
  140. #ifdef NO_LPC_FIX
  141.     /* Old decimation sometimes fails to recognise voiced. */
  142.     for (i = 0, j = 0; i < BUFLEN; i += DOWN)
  143.         d[j++] = w[i];
  144. #else
  145.     /* New: average rather than decimating. */
  146.     for (i = 0, j = 0; i < BUFLEN; j++) {
  147.         d[j] = 0;
  148.         for (rpos = 0; rpos < DOWN; rpos++) {
  149.         d[j] += w[i++];
  150.         }
  151.         /* d[j] /- DOWN;       Not actually necessary. */       
  152.     }
  153. #endif
  154.     auto_correl(d, BUFLEN / DOWN, PITCHORDER, r);
  155.     durbin(r, PITCHORDER, k, &g);
  156.     inverse_filter(d, k);
  157. #ifdef NO_LPC_FIX
  158.     auto_correl(d, BUFLEN / DOWN, MAXPER + 1, r);
  159. #else
  160.     auto_correl(d, BUFLEN / DOWN, MAXPER, r);
  161. #endif
  162.     rpos = 0;
  163.     rmax = 0.0;
  164.     for (i = MINPER; i <= MAXPER; i++) {
  165.         if (r[i] > rmax) {
  166.             rmax = r[i];
  167.             rpos = i;
  168.         }
  169.     }
  170.  
  171.     rm = r[rpos - 1];
  172.     rp = r[rpos + 1];
  173. #ifdef OLDWAY
  174.     rval = rmax / r[0];
  175. #endif
  176.  
  177.     a = 0.5 * rm - rmax + 0.5 * rp;
  178.     b = -0.5 * rm * (2.0 * rpos + 1.0) + 
  179.         2.0 * rpos * rmax + 0.5 * rp * (1.0 - 2.0 * rpos);
  180.     c = 0.5 * rm * (rpos * rpos + rpos) +
  181.         rmax * (1.0 - rpos * rpos) + 0.5 * rp * (rpos * rpos - rpos);
  182.  
  183.     x = -b / (2.0 * a);
  184.     y = a * x * x + b * x + c;
  185.     x *= DOWN;
  186.  
  187.     rmax = y;
  188. #ifdef OLDWAY
  189.     rval = rmax / r[0];
  190. #else
  191.     if (r[0] == 0.0) {
  192.         rval = 1.0;
  193.     } else {
  194.         rval = rmax / r[0];
  195.     }
  196. #endif
  197. #ifdef NO_LPC_FIX
  198.     if (rval >= 0.4 || (vuv == 3 && rval >= 0.3)) {
  199. #else
  200.     if ((rval >= 0.4 || (vuv == 3 && rval >= 0.3)) && (x > 0)) {
  201. #endif
  202.         *per = x;
  203.         vuv = (vuv & 1) * 2 + 1;
  204.     } else {
  205.         *per = 0.0;
  206.         vuv = (vuv & 1) * 2;
  207.     }
  208. }
  209.  
  210. #define s   work->w_s
  211. #define y   work->w_y
  212. #define h   work->w_h
  213.  
  214. void
  215. lpc_init(state)
  216.   lpcstate_t* state;
  217. {
  218.     int    i;
  219.     float    r, v, w, wcT;
  220.  
  221. #ifdef HEXDUMP
  222.         /*  Let's make sure the compiler hasn't inserted
  223.         padding in the lpcparams_t structure which will
  224.         make it incompatible with other machines.  */
  225.  
  226.     {
  227.         lpcparams_t *p = 0;
  228.  
  229.         if ((&p->gain != ((unsigned char *) 2)) ||
  230.         (p->k != ((signed char *) 4))) {
  231.                 fprintf(stderr, "Alignment problem in lpcparams.h structure.\n");
  232.                 fprintf(stderr, "Add definitions or compiler options in lpc directory\n");
  233.                 fprintf(stderr, "to guarantee this structure is packed.\n");
  234.         }
  235.     }
  236. #endif
  237.  
  238.     for (i = 0; i < BUFLEN; i++) {
  239.         s[i] = 0.0;
  240.         h[i] = WSCALE * (0.54 - 0.46 *
  241.                cos(2 * M_PI * i / (BUFLEN - 1.0)));
  242.     }
  243.     wcT = 2 * M_PI * FC / FS;
  244.     r = 0.36891079 * wcT;
  245.     v = 0.18445539 * wcT;
  246.     w = 0.92307712 * wcT;
  247.     fa[1] = -exp(-r);
  248.     fa[2] = 1.0 + fa[1];
  249.     fa[3] = -2.0 * exp(-v) * cos(w);
  250.     fa[4] = exp(-2.0 * v);
  251.     fa[5] = 1.0 + fa[3] + fa[4];
  252.  
  253.     u1 = 0.0;
  254.     yp1 = 0.0;
  255.     yp2 = 0.0;
  256.  
  257.     state->Oldper = 0.0;
  258.     state->OldG = 0.0;
  259.     for (i = 0; i <= LPC_FILTORDER; i++) {
  260.         state->Oldk[i] = 0.0;
  261.         state->bp[i] = 0.0;
  262.     }
  263.     state->pitchctr = 0;
  264. }
  265.  
  266. void
  267. lpc_analyze(buf, params)
  268.   const unsigned char *buf;
  269.   lpcparams_t *params;
  270. {
  271.     int    i, j;
  272. #define w   work->w_w
  273.     float    r[LPC_FILTORDER + 1];
  274.     float    per, G, k[LPC_FILTORDER + 1];
  275.  
  276.     for (i = 0, j = BUFLEN - FRAMESIZE; i < FRAMESIZE; i++, j++) {
  277.         s[j] = (float) (GAIN_ADJUST * ((ulaw_to_lin(*buf++)) / 32768.));
  278.         u = fa[2] * s[j] - fa[1] * u1;
  279.         y[j] = fa[5] * u1 - fa[3] * yp1 - fa[4] * yp2;
  280.         u1 = u;
  281.         yp2 = yp1;
  282.         yp1 = y[j];
  283.     }
  284.  
  285.     calc_pitch(y, &per);
  286.  
  287.     for (i = 0; i < BUFLEN; i++)
  288.         w[i] = s[i] * h[i];
  289.     auto_correl(w, BUFLEN, LPC_FILTORDER, r);
  290.     durbin(r, LPC_FILTORDER, k, &G);
  291.  
  292.     params->period = (unsigned short) htons((unsigned short) (per * (1 << 8)));
  293. #ifdef NO_LPC_FIX
  294.     params->gain = G * (1 << 8);
  295. #else
  296.     i = (int) (G * (1 << 8));
  297.     if (i > 255) {
  298.         i = 255;
  299.     }
  300.     params->gain = i;
  301. #endif
  302.     for (i = 0; i < LPC_FILTORDER; i++) {
  303. #ifdef NO_LPC_FIX
  304.         params->k[i] = k[i + 1] * (1 << 7);
  305. #else
  306.         float u = k[i + 1];
  307.  
  308.         if (u < -0.9999) {
  309.         u = -0.9999f;
  310.         } else if (u > 0.9999) {
  311.         u = 0.9999f;
  312.         }
  313.         params->k[i] = (signed char) (0.5 + (127.0 * k[i + 1]));
  314. #endif
  315.     }
  316.  
  317.     bcopy(s + FRAMESIZE, s, (BUFLEN - FRAMESIZE) * sizeof(s[0]));
  318.     bcopy(y + FRAMESIZE, y, (BUFLEN - FRAMESIZE) * sizeof(y[0]));
  319. #undef w
  320. }
  321.  
  322. void
  323. lpc_synthesize(buf, params, state)
  324.   unsigned char *buf;
  325.   lpcparams_t *params;
  326.   lpcstate_t* state;
  327. {
  328.     int i, j;
  329.     register double u, f, per, G, NewG, Ginc, Newper, perinc;
  330.     double k[LPC_FILTORDER + 1], Newk[LPC_FILTORDER + 1],
  331.            kinc[LPC_FILTORDER + 1];
  332.  
  333.     per = (double) ((unsigned short) ntohs(params->period)) / 256.;
  334.     G = (double) params->gain / 256.;
  335.     k[0] = 0.0;
  336.     for (i = 0; i < LPC_FILTORDER; i++)
  337.         k[i + 1] = (double) (params->k[i]) / 128.;
  338.  
  339.     G /= sqrt(BUFLEN / (per == 0.0? 3.0 : per));
  340.     Newper = state->Oldper;
  341.     NewG = state->OldG;
  342.     for (i = 1; i <= LPC_FILTORDER; i++)
  343.         Newk[i] = state->Oldk[i];
  344.  
  345.     if (state->Oldper != 0 && per != 0) {
  346.         perinc = (per - state->Oldper) / (double)FRAMESIZE;
  347.         Ginc = (G - state->OldG) / (double)FRAMESIZE;
  348.         for (i = 1; i <= LPC_FILTORDER; i++)
  349.             kinc[i] = (k[i] - state->Oldk[i]) / (double)FRAMESIZE;
  350.     } else {
  351.         perinc = 0.0;
  352.         Ginc = 0.0;
  353.         for (i = 1; i <= LPC_FILTORDER; i++)
  354.             kinc[i] = 0.0;
  355.     }
  356.  
  357.     if (Newper == 0)
  358.         state->pitchctr = 0;
  359.  
  360.     for (i = 0; i < FRAMESIZE; i++) {
  361.         if (Newper == 0) {
  362. #ifdef NO_LPC_FIX
  363.             u = ((double)random() / 2147483648.0) * NewG;
  364. #else
  365.             u = ((rand() / (1.0 + RAND_MAX)) - 0.5) * 1.5874 * NewG;
  366. #endif
  367.         } else {
  368.             if (state->pitchctr == 0) {
  369.                 u = NewG;
  370.                 state->pitchctr = (int) Newper;
  371.             } else {
  372.                 u = 0.0;
  373.                 state->pitchctr--;
  374.             }
  375.         }
  376.  
  377.         f = u;
  378.         for (j = LPC_FILTORDER; j >= 1; j--) {
  379.             register double b = state->bp[j - 1];
  380.             register double kj = Newk[j];
  381.             Newk[j] = kj + kinc[j];
  382.             f -= b * kj;
  383.             b += f * kj;
  384.             state->bp[j] = b;
  385.         }
  386.         state->bp[0] = f;
  387.  
  388. #ifdef NO_LPC_FIX
  389.         *buf++ = lin_to_ulaw((int)(f * 32768.0) & 0xffff);
  390. #else
  391.         u = f;
  392.         if (u < -0.9999) {
  393.             u = -0.9999;
  394.         } else if (u > 0.9999) {
  395.             u = 0.9999;
  396.         }
  397.         *buf++ = lin_to_ulaw((int) (u * 32767.0));
  398. #endif
  399.  
  400.         Newper += perinc;
  401.         NewG += Ginc;
  402.     }
  403.  
  404.     state->Oldper = per;
  405.     state->OldG = G;
  406.     for (i = 1; i <= LPC_FILTORDER; i++)
  407.         state->Oldk[i] = k[i];
  408. }
  409.  
  410. /*  LPC_START  --  Allocate working storage for LPC coder.  */
  411.  
  412. int lpc_start()
  413. {
  414.     if (work == NULL) {
  415.     work = (struct lpcwork *) malloc(sizeof(struct lpcwork));
  416.     }
  417.     return work != NULL;
  418. }
  419.  
  420. /*  LPC_END  --  Release working storage for LPC coder.  */
  421.  
  422. void lpc_end()
  423. {
  424.     if (work != NULL) {
  425.     free(work);
  426.     work = NULL;
  427.     }
  428. }
  429.